      subroutine NonResonantMatrix (o,Mbj2,dim_n,Cp,Cap,Rp,Rap,PE)
c
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer dim_tot,dim_een,dim_n,dim_PRC,dim_el
Cf2py intent(in) dim_n
c declaration of boundaries for the matrix
      integer c1_i,c1_e,c2_i,c2_e,c3_i,c3_e
c
c declaration of the matrix
c  output matrix
      real*8 o(1:dim_tot,1:dim_tot)
Cf2py intent(in,out) o
c  coefficient matrix (element wise squared)
      real*8 Mbj2(1:dim_een,1:dim_een)
Cf2py intent(in) Mbj2
c  capture part of the o matrix
      real*8 Ci(1:dim_een,1:dim_een),Cj(1:dim_een,1:dim_een)
      real*8 Cm11(1:dim_PRC*dim_el,1:dim_PRC*dim_el),Cm21(1:dim_PRC*dim_el,1:dim_PRC*dim_el)
      real*8 Cm12(1:dim_een,1:dim_een),Cm22(1:dim_een,1:dim_een)
c  recombination part of the o matrix
      real*8 Rj(1:dim_een,1:dim_een),Rk(1:dim_n, 1:dim_een)
c  photo excitation part of the o matrix
      real*8 PEi(1:dim_een,1:dim_een),PEk(1:dim_een,1:dim_een)
c  some basic matrix
      real*8 Sz(1:2,1:2),Sd(1:2,1:2),Id(1:dim_n,1:dim_n)
c
c declaration of some constant
      real*8 Cp,Cap,Calpha,Cbeta
Cf2py intent(in) Cp,Cap
      real*8 Rp,Rap
Cf2py intent(in) Rp,Rap
      real*8 PE
Cf2py intent(in) PE
c
c DEFINITIONS and INITIALIZATIONS
c constants
c  dimensions
      dim_PRC=2
      dim_el=2
      dim_tot=(dim_PRC*dim_el*dim_nu)+(dim_PRC*dim_el*dim_nu)+(dim_n)
      c1_i=1
      c1_e=(dim_PRC*dim_el*dim_n)
      c2_i=c1_e+1
      c2_e=c1_e+(dim_PRC*dim_el*dim_n)
      c3_i=c2_e+1
      c3_e=c2_e+dim_nu
c  base rates
      Calpha = (Cp - Cap) / 2.0
      Cbeta = (Cp + Cap) / 2.0
c matrix
      Ci=0
      Cj=0
      Cm11=0
      Cm21=0
      Cm12=0
      Cm22=0
      Ri=0
      Rj=0
      PEi=0
      PEk=0
      Sz=0
      Sd=0
      Id=0

c Construct self.NRrates
c Ci: assuming there is no coupling at this step (even if we write the transformation matrix)
      call kron(Cm11,Sz,2,2,Sz,2,2)
      call kron(Cm12,Cm11,4,4,Id,dim_n,dim_n)
      call kron(Cm21,Sd,2,2,Sd,2,2)
      call kron(Cm22,Cm21,4,4,Id,dim_n,dim_n)
      Ci = Calpha *  Cm12 + Cbeta * Cm22
        
      !Cj
      v=np.kron(np.array([Cp,Cap,Cap,Cp]),np.ones([dim_n]))
	do k=1,dim_een,1
      Cj(k,:)=transpose(matmul(v,Mbj2(:,k)))
	end do	

      !Rj
      R2=np.array([Rp,Rap,Rap,Rp])
      R3=np.ones([self.dim_n])
      R1=np.kron(R2,R3)
	do j=1,dim_een,1
      Rj(j,j)=matmul(transpose(R1),Mbj2(:,j))
	end do

      !Rk
	do j=1,dim_n,1
      R4=np.zeros([dim_n])
      R4(j)=1
      R4=np.kron(R2,R4)
	do i=1,dim_een,1
      Rk(j,i)=matmul(R4,Mbj2(:,i))
	end do
	end do
              
	!W: we assume there is no polarization of the excitation light ...
      PEi = PE/dim_een*np.kron(np.ones([1dim_PRC*dim_el,1]),Id)
      PEk = PE / dim_n * Id

      !Construct the rates matrix
      o(c1_i:c1_e, c1_i:c1_e) = -Ci
      o(c2_i:c2_e, c1_i:c1_e) = Cj
      o(c2_i:c2_e, c2_i:c2_e) = -Rj
      o(c3_i:c3_e, c2_i:c2_e) = Rk
      o(c3_i:c3_e, c3_i:c3_e) = -PEk
      o(c1_i:c1_e, c3_i:c3_e) = PEi

	end subroutine
